//
// File: timeFreqHRV.cpp
//
// MATLAB Coder version            : 5.4
// C/C++ source code generated on  : 23-Apr-2025 16:26:18
//

// Include Files
#include "timeFreqHRV.h"
#include "blockedSummation.h"
#include "colon.h"
#include "detrend.h"
#include "div.h"
#include "eml_mtimes_helper.h"
#include "fft.h"
#include "find.h"
#include "freqDomainHRV.h"
#include "freqDomainHRV_data.h"
#include "freqDomainHRV_initialize.h"
#include "freqDomainHRV_rtwutil.h"
#include "freqDomainHRV_types.h"
#include "hamming.h"
#include "ifft.h"
#include "interp1.h"
#include "linspace.h"
#include "lomb2.h"
#include "mean.h"
#include "minOrMax.h"
#include "pburg.h"
#include "rt_nonfinite.h"
#include "std.h"
#include "sum.h"
#include "trapz.h"
#include "var.h"
#include "wavelet.h"
#include "coder_array.h"
#include <cmath>
#include <string.h>

// Function Declarations
static void b_binary_expand_op(coder::array<double, 2U> &in1,
                               const struct9_T *in2);

static void binary_expand_op(struct9_T *in1,
                             const coder::array<double, 1U> &in2);

static void c_binary_expand_op(coder::array<double, 2U> &in1,
                               const struct9_T *in2);

static void calcAreas(const double F[449], const double PSD[449],
                      const double VLF[2], const double LF[2],
                      const double HF[2], struct2_T *output);

static void calcHRV(const coder::array<double, 1U> &F,
                    const coder::array<double, 2U> &PSD, const double VLF[2],
                    const double LF[2], const double HF[2], struct11_T *output);

static void calcLomb(const coder::array<double, 1U> &t,
                     const coder::array<double, 1U> &y, double nfft,
                     double maxF, double winSize, double overlap,
                     coder::array<double, 2U> &PSD, coder::array<double, 1U> &F,
                     coder::array<double, 1U> &T);

static double rt_powd_snf(double u0, double u1);

// Function Definitions
//
// Arguments    : coder::array<double, 2U> &in1
//                const struct9_T *in2
// Return Type  : void
//
static void b_binary_expand_op(coder::array<double, 2U> &in1,
                               const struct9_T *in2)
{
  int i;
  int loop_ub;
  int stride_0_1;
  int stride_1_1;
  if (in2->wav.hrv.aHF.size(1) == 1) {
    i = in2->wav.hrv.LFHF.size(1);
  } else {
    i = in2->wav.hrv.aHF.size(1);
  }
  in1.set_size(1, i);
  stride_0_1 = (in2->wav.hrv.LFHF.size(1) != 1);
  stride_1_1 = (in2->wav.hrv.aHF.size(1) != 1);
  if (in2->wav.hrv.aHF.size(1) == 1) {
    loop_ub = in2->wav.hrv.LFHF.size(1);
  } else {
    loop_ub = in2->wav.hrv.aHF.size(1);
  }
  for (i = 0; i < loop_ub; i++) {
    in1[i] =
        in2->wav.hrv.LFHF[i * stride_0_1] + in2->wav.hrv.aHF[i * stride_1_1];
  }
}

//
// Arguments    : struct9_T *in1
//                const coder::array<double, 1U> &in2
// Return Type  : void
//
static void binary_expand_op(struct9_T *in1,
                             const coder::array<double, 1U> &in2)
{
  coder::array<double, 1U> b_in1;
  int i;
  int loop_ub;
  int stride_0_0;
  int stride_1_0;
  if (in2.size(0) == 1) {
    i = in1->lomb.global.psd.size(0);
  } else {
    i = in2.size(0);
  }
  b_in1.set_size(i);
  stride_0_0 = (in1->lomb.global.psd.size(0) != 1);
  stride_1_0 = (in2.size(0) != 1);
  if (in2.size(0) == 1) {
    loop_ub = in1->lomb.global.psd.size(0);
  } else {
    loop_ub = in2.size(0);
  }
  for (i = 0; i < loop_ub; i++) {
    b_in1[i] = in1->lomb.global.psd[i * stride_0_0] * in2[i * stride_1_0];
  }
  in1->lomb.global.psd.set_size(b_in1.size(0));
  loop_ub = b_in1.size(0);
  for (i = 0; i < loop_ub; i++) {
    in1->lomb.global.psd[i] = b_in1[i];
  }
}

//
// Arguments    : coder::array<double, 2U> &in1
//                const struct9_T *in2
// Return Type  : void
//
static void c_binary_expand_op(coder::array<double, 2U> &in1,
                               const struct9_T *in2)
{
  int i;
  int loop_ub;
  int stride_0_1;
  int stride_1_1;
  int stride_2_1;
  if (in2->wav.hrv.aHF.size(1) == 1) {
    if (in2->wav.hrv.LFHF.size(1) == 1) {
      i = in2->wav.hrv.aVLF.size(1);
    } else {
      i = in2->wav.hrv.LFHF.size(1);
    }
  } else {
    i = in2->wav.hrv.aHF.size(1);
  }
  in1.set_size(1, i);
  stride_0_1 = (in2->wav.hrv.aVLF.size(1) != 1);
  stride_1_1 = (in2->wav.hrv.LFHF.size(1) != 1);
  stride_2_1 = (in2->wav.hrv.aHF.size(1) != 1);
  if (in2->wav.hrv.aHF.size(1) == 1) {
    if (in2->wav.hrv.LFHF.size(1) == 1) {
      loop_ub = in2->wav.hrv.aVLF.size(1);
    } else {
      loop_ub = in2->wav.hrv.LFHF.size(1);
    }
  } else {
    loop_ub = in2->wav.hrv.aHF.size(1);
  }
  for (i = 0; i < loop_ub; i++) {
    in1[i] = (in2->wav.hrv.aVLF[i * stride_0_1] +
              in2->wav.hrv.LFHF[i * stride_1_1]) +
             in2->wav.hrv.aHF[i * stride_2_1];
  }
}

//
// calcAreas - Calulates areas/energy under the PSD curve within the freq
// bands defined by VLF, LF, and HF. Returns areas/energies as ms^2,
// percentage, and normalized units. Also returns LF/HF ratio.
//
// Inputs:
//    PSD: PSD vector
//    F: Freq vector
//    VLF, LF, HF: array containing VLF, LF, and HF freq limits
//    flagNormalize: option to normalize PSD to max(PSD)
// Output:
//
// Usage:
//
//
//  Ref: This code is based on the calc_lfhf.m function from Gary Clifford's ECG
//  Toolbox.
//
// Arguments    : const double F[449]
//                const double PSD[449]
//                const double VLF[2]
//                const double LF[2]
//                const double HF[2]
//                struct2_T *output
// Return Type  : void
//
static void calcAreas(const double F[449], const double PSD[449],
                      const double VLF[2], const double LF[2],
                      const double HF[2], struct2_T *output)
{
  coder::array<double, 1U> b_PSD_data;
  coder::array<double, 1U> c_PSD_data;
  coder::array<double, 1U> d_PSD_data;
  double F_data[449];
  double PSD_data[449];
  double tmpF_data[449];
  double aHF;
  double aLF;
  double aTotal;
  double aVLF;
  double peakLF;
  double peakVLF;
  int F_size[2];
  int b_partialTrueCount;
  int iindx;
  int partialTrueCount;
  int trueCount;
  short b_tmp_data[449];
  short c_tmp_data[449];
  short d_tmp_data[449];
  short e_tmp_data[449];
  short f_tmp_data[449];
  short tmp_data[449];
  boolean_T bv[449];
  boolean_T bv1[449];
  boolean_T bv2[449];
  boolean_T bv3[449];
  boolean_T bv4[449];
  boolean_T bv5[449];
  // normalize PSD if needed
  //  find the indexes corresponding to the VLF, LF, and HF bands
  // Find peaks
  // VLF Peak
  partialTrueCount = 0;
  trueCount = 0;
  b_partialTrueCount = 0;
  aVLF = VLF[0];
  aTotal = VLF[1];
  aLF = LF[0];
  aHF = LF[1];
  peakVLF = HF[0];
  peakLF = HF[1];
  for (int i{0}; i < 449; i++) {
    double d;
    boolean_T b;
    boolean_T b1;
    d = F[i];
    b = (d >= aVLF);
    bv[i] = b;
    b1 = (d <= aTotal);
    bv1[i] = b1;
    bv2[i] = (d >= aLF);
    bv3[i] = (d <= aHF);
    bv4[i] = (d >= peakVLF);
    bv5[i] = (d <= peakLF);
    if (b && b1) {
      tmpF_data[partialTrueCount] = d;
      partialTrueCount++;
      trueCount++;
      tmp_data[b_partialTrueCount] = static_cast<short>(i + 1);
      b_partialTrueCount++;
    }
  }
  for (b_partialTrueCount = 0; b_partialTrueCount < trueCount;
       b_partialTrueCount++) {
    PSD_data[b_partialTrueCount] = PSD[tmp_data[b_partialTrueCount] - 1];
  }
  b_PSD_data.set(&PSD_data[0], trueCount);
  coder::internal::maximum(b_PSD_data, &aVLF, &iindx);
  peakVLF = tmpF_data[iindx - 1];
  // LF Peak
  partialTrueCount = 0;
  trueCount = 0;
  b_partialTrueCount = 0;
  for (int i{0}; i < 449; i++) {
    if (bv2[i] && bv3[i]) {
      tmpF_data[partialTrueCount] = F[i];
      partialTrueCount++;
      trueCount++;
      b_tmp_data[b_partialTrueCount] = static_cast<short>(i + 1);
      b_partialTrueCount++;
    }
  }
  for (b_partialTrueCount = 0; b_partialTrueCount < trueCount;
       b_partialTrueCount++) {
    PSD_data[b_partialTrueCount] = PSD[b_tmp_data[b_partialTrueCount] - 1];
  }
  c_PSD_data.set(&PSD_data[0], trueCount);
  coder::internal::maximum(c_PSD_data, &aVLF, &iindx);
  peakLF = tmpF_data[iindx - 1];
  // HF Peak
  partialTrueCount = 0;
  trueCount = 0;
  b_partialTrueCount = 0;
  for (int i{0}; i < 449; i++) {
    if (bv4[i] && bv5[i]) {
      tmpF_data[partialTrueCount] = F[i];
      partialTrueCount++;
      trueCount++;
      c_tmp_data[b_partialTrueCount] = static_cast<short>(i + 1);
      b_partialTrueCount++;
    }
  }
  for (b_partialTrueCount = 0; b_partialTrueCount < trueCount;
       b_partialTrueCount++) {
    PSD_data[b_partialTrueCount] = PSD[c_tmp_data[b_partialTrueCount] - 1];
  }
  d_PSD_data.set(&PSD_data[0], trueCount);
  coder::internal::maximum(d_PSD_data, &aVLF, &iindx);
  //  calculate raw areas (power under curve), within the freq bands (ms^2)
  trueCount = 0;
  partialTrueCount = 0;
  for (int i{0}; i < 449; i++) {
    if (bv[i] && bv1[i]) {
      trueCount++;
      d_tmp_data[partialTrueCount] = static_cast<short>(i + 1);
      partialTrueCount++;
    }
  }
  F_size[0] = 1;
  F_size[1] = trueCount;
  for (b_partialTrueCount = 0; b_partialTrueCount < trueCount;
       b_partialTrueCount++) {
    F_data[b_partialTrueCount] = F[d_tmp_data[b_partialTrueCount] - 1];
  }
  for (b_partialTrueCount = 0; b_partialTrueCount < trueCount;
       b_partialTrueCount++) {
    PSD_data[b_partialTrueCount] = PSD[d_tmp_data[b_partialTrueCount] - 1];
  }
  aVLF = coder::trapz(F_data, F_size, PSD_data, trueCount);
  trueCount = 0;
  partialTrueCount = 0;
  for (int i{0}; i < 449; i++) {
    if (bv2[i] && bv3[i]) {
      trueCount++;
      e_tmp_data[partialTrueCount] = static_cast<short>(i + 1);
      partialTrueCount++;
    }
  }
  F_size[0] = 1;
  F_size[1] = trueCount;
  for (b_partialTrueCount = 0; b_partialTrueCount < trueCount;
       b_partialTrueCount++) {
    F_data[b_partialTrueCount] = F[e_tmp_data[b_partialTrueCount] - 1];
  }
  for (b_partialTrueCount = 0; b_partialTrueCount < trueCount;
       b_partialTrueCount++) {
    PSD_data[b_partialTrueCount] = PSD[e_tmp_data[b_partialTrueCount] - 1];
  }
  aLF = coder::trapz(F_data, F_size, PSD_data, trueCount);
  trueCount = 0;
  partialTrueCount = 0;
  for (int i{0}; i < 449; i++) {
    if (bv4[i] && bv5[i]) {
      trueCount++;
      f_tmp_data[partialTrueCount] = static_cast<short>(i + 1);
      partialTrueCount++;
    }
  }
  F_size[0] = 1;
  F_size[1] = trueCount;
  for (b_partialTrueCount = 0; b_partialTrueCount < trueCount;
       b_partialTrueCount++) {
    F_data[b_partialTrueCount] = F[f_tmp_data[b_partialTrueCount] - 1];
  }
  for (b_partialTrueCount = 0; b_partialTrueCount < trueCount;
       b_partialTrueCount++) {
    PSD_data[b_partialTrueCount] = PSD[f_tmp_data[b_partialTrueCount] - 1];
  }
  aHF = coder::trapz(F_data, F_size, PSD_data, trueCount);
  // aVLF=aVLF/1e12;
  // aLF=aLF/1e12;
  // aHF=aHF/1e12;
  aTotal = (aVLF + aLF) + aHF;
  // calculate areas relative to the total area (%)
  // calculate normalized areas (relative to HF+LF, n.u.)
  // calculate LF/HF ratio
  // create output structure
  output->aVLF = std::round(aVLF * 100.0) / 100.0;
  //  round
  output->aLF = std::round(aLF * 100.0) / 100.0;
  output->aHF = std::round(aHF * 100.0) / 100.0;
  output->aTotal = std::round(aTotal * 100.0) / 100.0;
  output->pVLF = std::round(aVLF / aTotal * 100.0 * 10.0) / 10.0;
  output->pLF = std::round(aLF / aTotal * 100.0 * 10.0) / 10.0;
  output->pHF = std::round(aHF / aTotal * 100.0 * 10.0) / 10.0;
  aVLF = aLF + aHF;
  output->nLF = std::round(aLF / aVLF * 1000.0) / 1000.0;
  output->nHF = std::round(aHF / aVLF * 1000.0) / 1000.0;
  output->LFHF = std::round(aLF / aHF * 1000.0) / 1000.0;
  output->peakVLF = std::round(peakVLF * 100.0) / 100.0;
  output->peakLF = std::round(peakLF * 100.0) / 100.0;
  output->peakHF = std::round(tmpF_data[iindx - 1] * 100.0) / 100.0;
}

//
// calcAreas - Calulates areas/energy under the PSD curve within the freq
//  bands defined by VLF, LF, and HF. Returns areas/energies as ms^2,
//  percentage, and normalized units. Also returns LF/HF ratio.
//
//  Inputs:
//    PSD: PSD vector
//    F: Freq vector
//    VLF, LF, HF: array containing VLF, LF, and HF freq limits
//    flagVLF: flag to decide whether to calculate VLF hrv
//  Output:
//
//  Usage:
//
//
//  Ref: Modified from Gary Clifford's ECG Toolbox: calc_lfhf.m
//
// Arguments    : const coder::array<double, 1U> &F
//                const coder::array<double, 2U> &PSD
//                const double VLF[2]
//                const double LF[2]
//                const double HF[2]
//                struct11_T *output
// Return Type  : void
//
static void calcHRV(const coder::array<double, 1U> &F,
                    const coder::array<double, 2U> &PSD, const double VLF[2],
                    const double LF[2], const double HF[2], struct11_T *output)
{
  coder::array<double, 1U> b_PSD;
  coder::array<boolean_T, 1U> x;
  struct2_T expl_temp;
  int b_nz;
  int k;
  int nz;
  int vlen;
  output->aVLF.set_size(PSD.size(1));
  vlen = PSD.size(1);
  for (k = 0; k < vlen; k++) {
    output->aVLF[k] = 0.0;
  }
  output->aLF.set_size(PSD.size(1));
  vlen = PSD.size(1);
  for (k = 0; k < vlen; k++) {
    output->aLF[k] = 0.0;
  }
  output->aHF.set_size(PSD.size(1));
  vlen = PSD.size(1);
  for (k = 0; k < vlen; k++) {
    output->aHF[k] = 0.0;
  }
  output->aTotal.set_size(PSD.size(1));
  vlen = PSD.size(1);
  for (k = 0; k < vlen; k++) {
    output->aTotal[k] = 0.0;
  }
  output->pVLF.set_size(PSD.size(1));
  vlen = PSD.size(1);
  for (k = 0; k < vlen; k++) {
    output->pVLF[k] = 0.0;
  }
  output->pLF.set_size(PSD.size(1));
  vlen = PSD.size(1);
  for (k = 0; k < vlen; k++) {
    output->pLF[k] = 0.0;
  }
  output->pHF.set_size(PSD.size(1));
  vlen = PSD.size(1);
  for (k = 0; k < vlen; k++) {
    output->pHF[k] = 0.0;
  }
  output->nLF.set_size(PSD.size(1));
  vlen = PSD.size(1);
  for (k = 0; k < vlen; k++) {
    output->nLF[k] = 0.0;
  }
  output->nHF.set_size(PSD.size(1));
  vlen = PSD.size(1);
  for (k = 0; k < vlen; k++) {
    output->nHF[k] = 0.0;
  }
  output->LFHF.set_size(PSD.size(1));
  vlen = PSD.size(1);
  for (k = 0; k < vlen; k++) {
    output->LFHF[k] = 0.0;
  }
  output->peakVLF.set_size(PSD.size(1));
  vlen = PSD.size(1);
  for (k = 0; k < vlen; k++) {
    output->peakVLF[k] = 0.0;
  }
  output->peakLF.set_size(PSD.size(1));
  vlen = PSD.size(1);
  for (k = 0; k < vlen; k++) {
    output->peakLF[k] = 0.0;
  }
  output->peakHF.set_size(PSD.size(1));
  vlen = PSD.size(1);
  for (k = 0; k < vlen; k++) {
    output->peakHF[k] = 0.0;
  }
  k = PSD.size(1);
  for (nz = 0; nz < k; nz++) {
    vlen = PSD.size(0);
    b_PSD.set_size(PSD.size(0));
    for (b_nz = 0; b_nz < vlen; b_nz++) {
      b_PSD[b_nz] = PSD[b_nz + PSD.size(0) * nz];
    }
    calcAreas(F, b_PSD, VLF, LF, HF, &expl_temp);
    // create output structure
    output->aVLF[nz] = expl_temp.aVLF;
    output->aLF[nz] = expl_temp.aLF;
    output->aHF[nz] = expl_temp.aHF;
    output->aTotal[nz] = expl_temp.aTotal;
    output->pVLF[nz] = expl_temp.pVLF;
    output->pLF[nz] = expl_temp.pLF;
    output->pHF[nz] = expl_temp.pHF;
    output->nLF[nz] = expl_temp.nLF;
    output->nHF[nz] = expl_temp.nHF;
    output->LFHF[nz] = expl_temp.LFHF;
    output->peakVLF[nz] = expl_temp.peakVLF;
    output->peakLF[nz] = expl_temp.peakLF;
    output->peakHF[nz] = expl_temp.peakHF;
  }
  x.set_size(output->LFHF.size(0));
  vlen = output->LFHF.size(0);
  for (k = 0; k < vlen; k++) {
    x[k] = (output->LFHF[k] > 1.0);
  }
  vlen = x.size(0);
  if (x.size(0) == 0) {
    b_nz = 0;
  } else {
    b_nz = x[0];
    for (k = 2; k <= vlen; k++) {
      b_nz += x[k - 1];
    }
  }
  x.set_size(output->LFHF.size(0));
  vlen = output->LFHF.size(0);
  for (k = 0; k < vlen; k++) {
    x[k] = (output->LFHF[k] <= 1.0);
  }
  vlen = x.size(0);
  if (x.size(0) == 0) {
    nz = 0;
  } else {
    nz = x[0];
    for (k = 2; k <= vlen; k++) {
      nz += x[k - 1];
    }
  }
  output->rLFHF = static_cast<double>(b_nz) / static_cast<double>(nz);
}

//
// calLomb - Calculates PSD using windowed Lomb-Scargle method.
//
// Inputs:
// Outputs:
//
// Arguments    : const coder::array<double, 1U> &t
//                const coder::array<double, 1U> &y
//                double nfft
//                double maxF
//                double winSize
//                double overlap
//                coder::array<double, 2U> &PSD
//                coder::array<double, 1U> &F
//                coder::array<double, 1U> &T
// Return Type  : void
//
static void calcLomb(const coder::array<double, 1U> &t,
                     const coder::array<double, 1U> &y, double nfft,
                     double maxF, double winSize, double overlap,
                     coder::array<double, 2U> &PSD, coder::array<double, 1U> &F,
                     coder::array<double, 1U> &T)
{
  coder::array<double, 2U> iStart;
  coder::array<double, 2U> idx;
  coder::array<double, 1U> b_x;
  coder::array<double, 1U> r;
  coder::array<double, 1U> twt;
  coder::array<double, 1U> x;
  coder::array<double, 1U> y2;
  coder::array<double, 1U> z;
  double tau;
  double w;
  int b_loop_ub;
  int i;
  int i1;
  int loop_ub;
  unsigned int unnamed_idx_0;
  // get limits of windows
  if (t[t.size(0) - 1] >= winSize) {
    // windowArray - Seperates the input array s into smaller segments/windows
    // with lengths of winSize. winSize can be either in units of time or
    // samples. If winSize is given in units of time the first column of the
    // input array must containg time and method should be set to a value of 1.
    //
    // INPUTS:
    //    s: input array OR length of array. If outputing segs or using the time
    //        methods, the 1st dim must contain time values.
    //    winSize: size of segments to return
    //    method: method of segmenting array (0=by samples, 1=by time)
    //    overlap: amount that segments/windows overlap.
    //    flagCenters: decides whether to output center of windows instead of
    //    beginning
    //
    // OUTPUTS:
    //    idx: 2dim array that contains the start and end of each segment/window
    //    centers: array containg locations of center of windows (useful when
    //    using the time method)
    //
    // EXAMPLES:
    //    1. Seperate input using window size of 100 seconds and overlap
    //    of 50 seconds.
    //    idx=slidingWindow(input,100,50,1);
    //
    //    2. Seperate input using window size of 128 samples and overlap
    //    of 64 samples.
    //    idx=slidingWindow(input,128,64,0);
    //
    //    3. Seperate input using window size of 128 samples and overlap
    //    of 64 samples. Also outputs a cell array containing the actual
    //    segments [idx, segs]=slidingWindow(input,128,64,0,1);
    // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    //  Copyright (C) 2010, John T. Ramshur, jramshur@gmail.com
    //
    //  This file is part of HRVAS
    //
    //  HRVAS is free software: you can redistribute it and/or modify
    //  it under the terms of the GNU General Public License as published by
    //  the Free Software Foundation, either version 3 of the License, or
    //  (at your option) any later version.
    //
    //  HRVAS is distributed in the hope that it will be useful,
    //  but WITHOUT ANY WARRANTY; without even the implied warranty of
    //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    //  GNU General Public License for more details.
    //
    //  You should have received a copy of the GNU General Public License
    //  along with Foobar.  If not, see <http://www.gnu.org/licenses/>.
    // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    // Check inputs
    //  if nargin<5
    //      flagCenters=false;
    //      c=[];
    //      %includeSeg=0;
    //  elseif nargin<4
    //      overlap=0;
    //      includeSeg=0;
    //  elseif nargin<3
    //      overlap=0;
    //      method=0;
    //      includeSeg=0;
    //  elseif nargin<3
    //      error('slidingWindow: Too few input arguments.')
    //      return;
    //  end
    // check inputs
    //  if overlap>=winSize
    //      error('slidingWindow: overlap cannot be >= winSize.')
    //  end
    //  [m,n]=size(s);
    //  if (m==1) && (n==1) && (method==1)
    //      error('slidingWindow: Invalid input array for this windowing
    //      method.') return;
    //  end
    tau = winSize - overlap;
    w = (static_cast<double>(t.size(0)) - winSize) + 1.0;
    if (std::isnan(tau)) {
      iStart.set_size(1, 1);
      iStart[0] = rtNaN;
    } else if ((tau == 0.0) || ((w > 1.0) && (tau < 0.0)) ||
               ((w < 1.0) && (tau > 0.0))) {
      iStart.set_size(1, 0);
    } else if (std::isinf(w) && std::isinf(tau)) {
      iStart.set_size(1, 1);
      iStart[0] = rtNaN;
    } else if (std::isinf(tau)) {
      iStart.set_size(1, 1);
      iStart[0] = 1.0;
    } else if (std::floor(tau) == tau) {
      loop_ub = static_cast<int>((w - 1.0) / tau);
      iStart.set_size(1, loop_ub + 1);
      for (i = 0; i <= loop_ub; i++) {
        iStart[i] = tau * static_cast<double>(i) + 1.0;
      }
    } else {
      coder::eml_float_colon(tau, w, iStart);
    }
    //      nseg = length(iStart);
    idx.set_size(iStart.size(1), 2);
    loop_ub = iStart.size(1);
    for (i = 0; i < loop_ub; i++) {
      idx[i] = iStart[i];
    }
    loop_ub = iStart.size(1);
    for (i = 0; i < loop_ub; i++) {
      idx[i + idx.size(0)] = iStart[i] + (winSize - 1.0);
    }
    //  if includeSeg
    //      segs=cell(size(idx,1));
    //      for i=1:size(idx,1)
    //          segs{i}=s(idx(i,1):idx(i,2));%code here for outputing
    //      end
    //  end
  } else {
    idx.set_size(1, 2);
    idx[0] = 1.0;
    idx[1] = t.size(0);
  }
  tau = std::round(winSize / 2.0);
  loop_ub = idx.size(0);
  T.set_size(idx.size(0));
  for (i = 0; i < loop_ub; i++) {
    T[i] = t[static_cast<int>(idx[i]) - 1] + tau;
  }
  // estimate the center of the windows for plotting
  // preallocate memory
  // number of PSD/windows
  PSD.set_size(static_cast<int>(nfft), idx.size(0));
  coder::linspace(maxF - maxF / nfft, nfft, iStart);
  F.set_size(iStart.size(1));
  loop_ub = iStart.size(1);
  for (i = 0; i < loop_ub; i++) {
    F[i] = iStart[i];
  }
  i = idx.size(0);
  if (idx.size(0) - 1 >= 0) {
    unnamed_idx_0 = static_cast<unsigned int>(F.size(0));
    b_loop_ub = F.size(0);
    i1 = F.size(0);
  }
  for (int b_i{0}; b_i < i; b_i++) {
    int i2;
    int i3;
    int k;
    int nx;
    // Prepare y2 and t2
    tau = idx[b_i];
    w = idx[b_i + idx.size(0)];
    if (tau > w) {
      k = 0;
      nx = 0;
      i2 = 0;
      i3 = 0;
    } else {
      k = static_cast<int>(tau) - 1;
      nx = static_cast<int>(w);
      i2 = static_cast<int>(tau) - 1;
      i3 = static_cast<int>(w);
    }
    // remove linear trend
    loop_ub = nx - k;
    y2.set_size(loop_ub);
    for (nx = 0; nx < loop_ub; nx++) {
      y2[nx] = y[k + nx];
    }
    coder::detrend(y2);
    tau = coder::blockedSummation(y2, y2.size(0)) /
          static_cast<double>(y2.size(0));
    loop_ub = y2.size(0);
    for (k = 0; k < loop_ub; k++) {
      y2[k] = y2[k] - tau;
    }
    // remove mean
    // Calculate un-normalized lomb PSD
    //
    //   [Pn] = lomb(t, y, f)
    //
    //   Uses Lomb's method to compute normalized
    //   periodogram values "Pn" as a function of
    //   supplied vector of frequencies "f" for
    //   input vectors "t" (time) and "y" (observations).
    //   Also returned is probability "Prob" of same
    //   length as Pn (and f) that the null hypothesis
    //   is valid. If f is not supplied it assumes
    //   f =  [1/1024 : 1/1024 : 0.5/min(diff(t))];
    //   x and y must be the same length.
    //  See also:
    //
    //  [1] N.R. Lomb, ``Least-squares frequency analysis of
    //  unequally spaced data,'' Astrophysical and Space Science,
    //  (39) pp. 447--462, 1976.   ... and
    //
    //  [2] J.~D. Scargle, ``Studies in astronomical time series analysis.
    //  II. Statistical aspects of spectral analysis of unevenly spaced data,''
    //  Astrophysical Journal, vol. 263, pp. 835--853, 1982.
    //
    //  [3] T. Thong, "Lomb-Welch Periodogram for Non-uniform Sampling",
    //  Proceedings for the 26th anual international conference of the IEEE
    //  EMBS, Sept 1-5, 2004. check inputs subtract mean, compute variance,
    //  initialize Pn
    tau = coder::blockedSummation(y2, y2.size(0)) /
          static_cast<double>(y2.size(0));
    z.set_size(y2.size(0));
    loop_ub = y2.size(0);
    for (k = 0; k < loop_ub; k++) {
      z[k] = y2[k] - tau;
    }
    r.set_size(static_cast<int>(unnamed_idx_0));
    for (k = 0; k < b_loop_ub; k++) {
      r[k] = 0.0;
    }
    // 	now do main loop for all frequencies
    for (int c_i{0}; c_i < i1; c_i++) {
      w = 6.2831853071795862 * F[c_i];
      if (w > 0.0) {
        double a;
        a = 2.0 * w;
        loop_ub = i3 - i2;
        twt.set_size(loop_ub);
        for (k = 0; k < loop_ub; k++) {
          twt[k] = a * t[i2 + k];
        }
        x.set_size(twt.size(0));
        nx = twt.size(0);
        for (k = 0; k < nx; k++) {
          x[k] = twt[k];
        }
        nx = twt.size(0);
        for (k = 0; k < nx; k++) {
          x[k] = std::sin(x[k]);
        }
        nx = twt.size(0);
        for (k = 0; k < nx; k++) {
          twt[k] = std::cos(twt[k]);
        }
        tau = rt_atan2d_snf(coder::blockedSummation(x, x.size(0)),
                            coder::blockedSummation(twt, twt.size(0))) /
              2.0 / w;
        twt.set_size(loop_ub);
        for (k = 0; k < loop_ub; k++) {
          twt[k] = w * (t[i2 + k] - tau);
        }
        x.set_size(twt.size(0));
        loop_ub = twt.size(0);
        for (k = 0; k < loop_ub; k++) {
          x[k] = twt[k];
        }
        nx = twt.size(0);
        for (k = 0; k < nx; k++) {
          x[k] = std::cos(x[k]);
        }
        b_x.set_size(twt.size(0));
        loop_ub = twt.size(0);
        for (k = 0; k < loop_ub; k++) {
          b_x[k] = twt[k];
        }
        nx = twt.size(0);
        for (k = 0; k < nx; k++) {
          b_x[k] = std::sin(b_x[k]);
        }
        loop_ub = z.size(0);
        if (z.size(0) == x.size(0)) {
          x.set_size(z.size(0));
          for (k = 0; k < loop_ub; k++) {
            x[k] = z[k] * x[k];
          }
        } else {
          times(x, z);
        }
        a = coder::blockedSummation(x, x.size(0));
        loop_ub = z.size(0);
        if (z.size(0) == b_x.size(0)) {
          b_x.set_size(z.size(0));
          for (k = 0; k < loop_ub; k++) {
            b_x[k] = z[k] * b_x[k];
          }
        } else {
          times(b_x, z);
        }
        w = coder::blockedSummation(b_x, b_x.size(0));
        x.set_size(twt.size(0));
        loop_ub = twt.size(0);
        for (k = 0; k < loop_ub; k++) {
          x[k] = twt[k];
        }
        nx = twt.size(0);
        for (k = 0; k < nx; k++) {
          x[k] = std::cos(x[k]);
        }
        nx = twt.size(0);
        for (k = 0; k < nx; k++) {
          twt[k] = std::sin(twt[k]);
        }
        loop_ub = x.size(0);
        for (k = 0; k < loop_ub; k++) {
          tau = x[k];
          x[k] = tau * tau;
        }
        b_x.set_size(twt.size(0));
        loop_ub = twt.size(0);
        for (k = 0; k < loop_ub; k++) {
          tau = twt[k];
          b_x[k] = tau * tau;
        }
        r[c_i] = a * a / coder::blockedSummation(x, x.size(0)) +
                 w * w / coder::blockedSummation(b_x, b_x.size(0));
      } else {
        double a;
        loop_ub = z.size(0);
        nx = i3 - i2;
        if (z.size(0) == nx) {
          x.set_size(z.size(0));
          for (k = 0; k < loop_ub; k++) {
            x[k] = z[k] * t[i2 + k];
          }
        } else {
          binary_expand_op(x, z, t, i2, i3 - 1);
        }
        a = coder::blockedSummation(x, x.size(0));
        x.set_size(nx);
        for (k = 0; k < nx; k++) {
          tau = t[i2 + k];
          x[k] = tau * tau;
        }
        r[c_i] = a * a / coder::blockedSummation(x, x.size(0));
      }
    }
    //  return denormalized spectrum (see T. Thong)
    loop_ub = r.size(0);
    for (k = 0; k < loop_ub; k++) {
      PSD[k + PSD.size(0) * b_i] = r[k] / static_cast<double>(y2.size(0));
    }
  }
}

//
// Arguments    : double u0
//                double u1
// Return Type  : double
//
static double rt_powd_snf(double u0, double u1)
{
  double y;
  if (std::isnan(u0) || std::isnan(u1)) {
    y = rtNaN;
  } else {
    double d;
    double d1;
    d = std::abs(u0);
    d1 = std::abs(u1);
    if (std::isinf(u1)) {
      if (d == 1.0) {
        y = 1.0;
      } else if (d > 1.0) {
        if (u1 > 0.0) {
          y = rtInf;
        } else {
          y = 0.0;
        }
      } else if (u1 > 0.0) {
        y = 0.0;
      } else {
        y = rtInf;
      }
    } else if (d1 == 0.0) {
      y = 1.0;
    } else if (d1 == 1.0) {
      if (u1 > 0.0) {
        y = u0;
      } else {
        y = 1.0 / u0;
      }
    } else if (u1 == 2.0) {
      y = u0 * u0;
    } else if ((u1 == 0.5) && (u0 >= 0.0)) {
      y = std::sqrt(u0);
    } else if ((u0 < 0.0) && (u1 > std::floor(u1))) {
      y = rtNaN;
    } else {
      y = std::pow(u0, u1);
    }
  }
  return y;
}

//
// timeFreqHRV - calculates time-freq HRV using ar, lomb, and CWT methods
//
//  Inputs:   ibi = 2Dim array of [time (s) ,inter-beat interval (s)]
//            nibi = IBI with trend still present (non-detrended). Used for CWT.
//            VLF,LF,HF = arrays containing limits of VLF, LF, and HF freq.
//            bands winSize = # of samples in window noverlap = # of samples to
//            overlap fs = cubic spline interpolation rate / resample rate (Hz)
//            nfft = # of points in the frequency axis
//            methods = cell array containing UP TO three strings that tell
//            the function what methods to include in calculating HRV
//            {'ar','lomb','wavelet'}
//  Outputs:  output is a structure containg all HRV. One field for each PSD
//  method
//            Output units include:
//                peakHF,peakLF,peakVLF (Hz)
//                aHF,aLF,aVLF (ms^2)
//                pHF,pLF,pVLF (%)
//                nHF,nLF,nVLF (%)
//                lfhf,rlfhf
//                PSD (ms^2/Hz)
//                F (Hz)
//                T (s)
//  Usage:   n/a
//
// Arguments    : coder::array<double, 2U> &ibi
//                const double VLF[2]
//                const double LF[2]
//                const double HF[2]
//                double AR_order
//                double winSize
//                double overlap
//                double nfft
//                double fs
//                struct9_T *output
// Return Type  : void
//
void timeFreqHRV(coder::array<double, 2U> &ibi, const double VLF[2],
                 const double LF[2], const double HF[2], double AR_order,
                 double winSize, double overlap, double nfft, double fs,
                 struct9_T *output)
{
  static const double b[449]{1.0,
                             1.0108892860517005,
                             1.0218971486541166,
                             1.0330248790212284,
                             1.0442737824274138,
                             1.0556451783605572,
                             1.0671404006768237,
                             1.0787607977571199,
                             1.0905077326652577,
                             1.1023825833078409,
                             1.1143867425958924,
                             1.1265216186082418,
                             1.1387886347566916,
                             1.1511892299529827,
                             1.1637248587775775,
                             1.1763969916502812,
                             1.189207115002721,
                             1.2021567314527031,
                             1.215247359980469,
                             1.22848053610687,
                             1.241857812073484,
                             1.2553807570246911,
                             1.2690509571917332,
                             1.2828700160787783,
                             1.2968395546510096,
                             1.3109612115247644,
                             1.3252366431597413,
                             1.3396675240533029,
                             1.3542555469368927,
                             1.3690024229745905,
                             1.383909881963832,
                             1.398979672538311,
                             1.4142135623730951,
                             1.42961333839197,
                             1.4451808069770467,
                             1.4609177941806468,
                             1.4768261459394993,
                             1.4929077282912648,
                             1.5091644275934226,
                             1.5255981507445382,
                             1.542210825407941,
                             1.5590044002378369,
                             1.5759808451078865,
                             1.593142151342267,
                             1.6104903319492543,
                             1.6280274218573478,
                             1.6457554781539647,
                             1.6636765803267364,
                             1.6817928305074292,
                             1.7001063537185235,
                             1.7186192981224779,
                             1.7373338352737062,
                             1.7562521603732995,
                             1.7753764925265212,
                             1.7947090750031072,
                             1.8142521755003989,
                             1.8340080864093424,
                             1.8539791250833857,
                             1.8741676341103,
                             1.8945759815869656,
                             1.9152065613971474,
                             1.9360617934922943,
                             1.9571441241754002,
                             1.9784560263879509,
                             2.0,
                             2.021778572103401,
                             2.0437942973082333,
                             2.0660497580424568,
                             2.0885475648548275,
                             2.1112903567211143,
                             2.1342808013536474,
                             2.1575215955142397,
                             2.1810154653305154,
                             2.2047651666156818,
                             2.2287734851917849,
                             2.2530432372164837,
                             2.2775772695133831,
                             2.3023784599059653,
                             2.327449717555155,
                             2.3527939833005624,
                             2.3784142300054421,
                             2.4043134629054062,
                             2.4304947199609379,
                             2.45696107221374,
                             2.483715624146968,
                             2.5107615140493822,
                             2.5381019143834664,
                             2.5657400321575565,
                             2.5936791093020193,
                             2.6219224230495288,
                             2.6504732863194826,
                             2.6793350481066058,
                             2.7085110938737853,
                             2.738004845949181,
                             2.767819763927664,
                             2.797959345076622,
                             2.82842712474619,
                             2.85922667678394,
                             2.8903616139540933,
                             2.9218355883612936,
                             2.9536522918789987,
                             2.9858154565825297,
                             3.0183288551868452,
                             3.0511963014890764,
                             3.0844216508158819,
                             3.1180088004756739,
                             3.151961690215773,
                             3.186284302684534,
                             3.2209806638985086,
                             3.2560548437146957,
                             3.2915109563079294,
                             3.3273531606534728,
                             3.3635856610148585,
                             3.400212707437047,
                             3.4372385962449559,
                             3.4746676705474124,
                             3.5125043207465989,
                             3.5507529850530424,
                             3.5894181500062143,
                             3.6285043510007977,
                             3.6680161728186849,
                             3.7079582501667714,
                             3.7483352682206,
                             3.7891519631739312,
                             3.8304131227942948,
                             3.8721235869845887,
                             3.9142882483508004,
                             3.9569120527759019,
                             4.0,
                             4.0435571442068019,
                             4.0875885946164665,
                             4.1320995160849137,
                             4.177095129709655,
                             4.2225807134422286,
                             4.2685616027072948,
                             4.3150431910284794,
                             4.3620309306610308,
                             4.4095303332313636,
                             4.45754697038357,
                             4.5060864744329674,
                             4.5551545390267663,
                             4.6047569198119307,
                             4.65489943511031,
                             4.7055879666011249,
                             4.7568284600108841,
                             4.8086269258108123,
                             4.8609894399218758,
                             4.91392214442748,
                             4.967431248293936,
                             5.0215230280987644,
                             5.0762038287669329,
                             5.1314800643151131,
                             5.1873582186040386,
                             5.2438448460990577,
                             5.3009465726389653,
                             5.3586700962132117,
                             5.4170221877475706,
                             5.4760096918983621,
                             5.5356395278553281,
                             5.5959186901532441,
                             5.65685424949238,
                             5.71845335356788,
                             5.7807232279081866,
                             5.8436711767225873,
                             5.9073045837579974,
                             5.9716309131650593,
                             6.0366577103736905,
                             6.1023926029781528,
                             6.1688433016317639,
                             6.2360176009513477,
                             6.303923380431546,
                             6.372568605369068,
                             6.4419613277970171,
                             6.5121096874293913,
                             6.5830219126158589,
                             6.6547063213069455,
                             6.7271713220297169,
                             6.8004254148740939,
                             6.8744771924899117,
                             6.9493353410948249,
                             7.0250086414931978,
                             7.1015059701060848,
                             7.1788363000124287,
                             7.2570087020015954,
                             7.33603234563737,
                             7.4159165003335428,
                             7.4966705364412,
                             7.5783039263478624,
                             7.66082624558859,
                             7.7442471739691774,
                             7.8285764967016007,
                             7.9138241055518037,
                             8.0,
                             8.0871142884136038,
                             8.175177189232933,
                             8.2641990321698273,
                             8.35419025941931,
                             8.4451614268844573,
                             8.53712320541459,
                             8.6300863820569589,
                             8.7240618613220615,
                             8.8190606664627271,
                             8.91509394076714,
                             9.0121729488659348,
                             9.1103090780535325,
                             9.2095138396238614,
                             9.30979887022062,
                             9.41117593320225,
                             9.5136569200217682,
                             9.6172538516216246,
                             9.7219788798437516,
                             9.82784428885496,
                             9.934862496587872,
                             10.043046056197529,
                             10.152407657533866,
                             10.262960128630226,
                             10.374716437208077,
                             10.487689692198115,
                             10.601893145277931,
                             10.717340192426423,
                             10.834044375495141,
                             10.952019383796724,
                             11.071279055710656,
                             11.191837380306488,
                             11.313708498984759,
                             11.43690670713576,
                             11.561446455816373,
                             11.687342353445175,
                             11.814609167515995,
                             11.943261826330119,
                             12.073315420747381,
                             12.204785205956306,
                             12.337686603263528,
                             12.472035201902695,
                             12.607846760863092,
                             12.745137210738136,
                             12.883922655594034,
                             13.024219374858783,
                             13.166043825231718,
                             13.309412642613891,
                             13.454342644059434,
                             13.600850829748188,
                             13.748954384979823,
                             13.89867068218965,
                             14.050017282986396,
                             14.20301194021217,
                             14.357672600024857,
                             14.514017404003191,
                             14.672064691274739,
                             14.831833000667086,
                             14.9933410728824,
                             15.156607852695725,
                             15.321652491177179,
                             15.488494347938355,
                             15.657152993403201,
                             15.827648211103607,
                             16.0,
                             16.174228576827208,
                             16.350354378465866,
                             16.528398064339655,
                             16.70838051883862,
                             16.890322853768915,
                             17.074246410829179,
                             17.260172764113918,
                             17.448123722644123,
                             17.638121332925454,
                             17.830187881534279,
                             18.02434589773187,
                             18.220618156107065,
                             18.419027679247723,
                             18.61959774044124,
                             18.8223518664045,
                             19.027313840043536,
                             19.234507703243249,
                             19.443957759687503,
                             19.65568857770992,
                             19.869724993175744,
                             20.086092112395058,
                             20.304815315067732,
                             20.525920257260452,
                             20.749432874416154,
                             20.975379384396231,
                             21.203786290555861,
                             21.434680384852847,
                             21.668088750990282,
                             21.904038767593448,
                             22.142558111421312,
                             22.383674760612976,
                             22.627416997969519,
                             22.87381341427152,
                             23.122892911632746,
                             23.374684706890349,
                             23.62921833503199,
                             23.886523652660237,
                             24.146630841494762,
                             24.409570411912611,
                             24.675373206527055,
                             24.944070403805391,
                             25.215693521726184,
                             25.490274421476272,
                             25.767845311188069,
                             26.048438749717565,
                             26.332087650463436,
                             26.618825285227782,
                             26.908685288118868,
                             27.201701659496376,
                             27.497908769959647,
                             27.7973413643793,
                             28.100034565972791,
                             28.406023880424339,
                             28.715345200049715,
                             29.028034808006382,
                             29.344129382549479,
                             29.663666001334171,
                             29.9866821457648,
                             30.31321570539145,
                             30.643304982354358,
                             30.97698869587671,
                             31.314305986806403,
                             31.655296422207215,
                             32.0,
                             32.348457153654415,
                             32.700708756931732,
                             33.056796128679309,
                             33.41676103767724,
                             33.780645707537829,
                             34.148492821658358,
                             34.520345528227836,
                             34.896247445288246,
                             35.276242665850909,
                             35.660375763068558,
                             36.048691795463739,
                             36.44123631221413,
                             36.838055358495446,
                             37.239195480882479,
                             37.644703732809,
                             38.054627680087073,
                             38.4690154064865,
                             38.887915519375007,
                             39.311377155419841,
                             39.739449986351488,
                             40.172184224790115,
                             40.609630630135463,
                             41.051840514520904,
                             41.498865748832308,
                             41.950758768792461,
                             42.407572581111722,
                             42.869360769705693,
                             43.336177501980565,
                             43.8080775351869,
                             44.285116222842625,
                             44.767349521225952,
                             45.254833995939038,
                             45.747626828543041,
                             46.245785823265493,
                             46.7493694137807,
                             47.258436670063979,
                             47.773047305320475,
                             48.293261682989524,
                             48.819140823825222,
                             49.350746413054111,
                             49.888140807610782,
                             50.431387043452368,
                             50.980548842952544,
                             51.535690622376137,
                             52.096877499435131,
                             52.664175300926871,
                             53.237650570455564,
                             53.817370576237735,
                             54.403403318992751,
                             54.995817539919294,
                             55.5946827287586,
                             56.200069131945583,
                             56.812047760848678,
                             57.430690400099429,
                             58.056069616012763,
                             58.688258765098958,
                             59.327332002668342,
                             59.9733642915296,
                             60.6264314107829,
                             61.286609964708717,
                             61.953977391753419,
                             62.628611973612806,
                             63.31059284441443,
                             64.0,
                             64.69691430730883,
                             65.401417513863464,
                             66.113592257358619,
                             66.83352207535448,
                             67.561291415075658,
                             68.296985643316717,
                             69.040691056455671,
                             69.792494890576492,
                             70.552485331701817,
                             71.320751526137116,
                             72.097383590927478,
                             72.88247262442826,
                             73.676110716990891,
                             74.478390961764958,
                             75.289407465618,
                             76.109255360174146,
                             76.938030812973,
                             77.775831038750013,
                             78.622754310839682,
                             79.478899972702976,
                             80.34436844958023,
                             81.219261260270926,
                             82.103681029041809,
                             82.997731497664617,
                             83.901517537584922,
                             84.815145162223445,
                             85.738721539411387,
                             86.67235500396113,
                             87.6161550703738,
                             88.570232445685249,
                             89.5346990424519,
                             90.509667991878075,
                             91.495253657086081,
                             92.491571646530986,
                             93.4987388275614,
                             94.516873340127958,
                             95.546094610640949,
                             96.586523365979048,
                             97.638281647650444,
                             98.701492826108222,
                             99.776281615221563,
                             100.86277408690474,
                             101.96109768590509,
                             103.07138124475227,
                             104.19375499887026,
                             105.32835060185374,
                             106.47530114091113,
                             107.63474115247547,
                             108.8068066379855,
                             109.99163507983859,
                             111.1893654575172,
                             112.40013826389117,
                             113.62409552169736,
                             114.86138080019886,
                             116.11213923202553,
                             117.37651753019792,
                             118.65466400533668,
                             119.9467285830592,
                             121.2528628215658,
                             122.57321992941743,
                             123.90795478350684,
                             125.25722394722561,
                             126.62118568882886,
                             128.0};
  coder::array<creal_T, 2U> b_f;
  coder::array<creal_T, 2U> f;
  coder::array<creal_T, 2U> r2;
  coder::array<creal_T, 2U> wave;
  coder::array<double, 2U> b_scale_avg;
  coder::array<double, 2U> iStart;
  coder::array<double, 2U> idx;
  coder::array<double, 2U> r;
  coder::array<double, 2U> r1;
  coder::array<double, 2U> scale_avg;
  coder::array<double, 1U> b_ibi;
  coder::array<double, 1U> c_ibi;
  coder::array<double, 1U> psd;
  coder::array<boolean_T, 2U> x;
  double b_d;
  double d;
  double d_tmp;
  double variance;
  double xtmp;
  int tmp_data[449];
  int tmp_size[2];
  int a1;
  int i;
  int loop_ub;
  int nx;
  int y;
  boolean_T b_output[449];
  if (!isInitialized_freqDomainHRV) {
    freqDomainHRV_initialize();
  }
  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  //  Copyright (C) 2010, John T. Ramshur, jramshur@gmail.com
  //
  //  This file is part of HRVAS
  //
  //  HRVAS is free software: you can redistribute it and/or modify
  //  it under the terms of the GNU General Public License as published by
  //  the Free Software Foundation, either version 3 of the License, or
  //  (at your option) any later version.
  //
  //  HRVAS is distributed in the hope that it will be useful,
  //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  //  GNU General Public License for more details.
  //
  //  You should have received a copy of the GNU General Public License
  //  along with Foobar.  If not, see <http://www.gnu.org/licenses/>.
  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  // check input
  //      if nargin<10; error('Not enough input arguments!'); end
  //      if nargin<11; methods={'ar','lomb','wavelet'}; end
  //
  // assumes ibi units are seconds
  nx = ibi.size(0) - 1;
  b_ibi.set_size(ibi.size(0));
  for (i = 0; i <= nx; i++) {
    b_ibi[i] = ibi[i + ibi.size(0)] * 1000.0;
  }
  loop_ub = b_ibi.size(0);
  for (i = 0; i < loop_ub; i++) {
    ibi[i + ibi.size(0)] = b_ibi[i];
  }
  // convert ibi units from s to ms
  //      nibi(:,2)=nibi(:,2).*1000; %convert ibi units from s to ms
  // time
  //      clear ibi; %don't need it anymore
  // AR
  //      output.ar.t=t;
  // calAR - Calculates PSD using windowed Burg method.
  //
  // Inputs:
  // Outputs:
  xtmp = winSize * fs;
  //  (samples)
  // resample
  d_tmp = 1.0 / fs;
  if (std::isnan(ibi[0]) || std::isnan(d_tmp) ||
      std::isnan(ibi[ibi.size(0) - 1])) {
    output->wav.hrv.aVLF.set_size(1, 1);
    output->wav.hrv.aVLF[0] = rtNaN;
  } else if ((d_tmp == 0.0) ||
             ((ibi[0] < ibi[ibi.size(0) - 1]) && (d_tmp < 0.0)) ||
             ((ibi[ibi.size(0) - 1] < ibi[0]) && (d_tmp > 0.0))) {
    output->wav.hrv.aVLF.set_size(1, 0);
  } else if ((std::isinf(ibi[0]) || std::isinf(ibi[ibi.size(0) - 1])) &&
             (std::isinf(d_tmp) || (ibi[0] == ibi[ibi.size(0) - 1]))) {
    output->wav.hrv.aVLF.set_size(1, 1);
    output->wav.hrv.aVLF[0] = rtNaN;
  } else if (std::isinf(d_tmp)) {
    output->wav.hrv.aVLF.set_size(1, 1);
    output->wav.hrv.aVLF[0] = ibi[0];
  } else if ((std::floor(ibi[0]) == ibi[0]) && (std::floor(d_tmp) == d_tmp)) {
    output->wav.hrv.aVLF.set_size(
        1, static_cast<int>((ibi[ibi.size(0) - 1] - ibi[0]) / d_tmp) + 1);
    loop_ub = static_cast<int>((ibi[ibi.size(0) - 1] - ibi[0]) / d_tmp);
    for (i = 0; i <= loop_ub; i++) {
      output->wav.hrv.aVLF[i] = ibi[0] + d_tmp * static_cast<double>(i);
    }
  } else {
    coder::eml_float_colon(ibi[0], d_tmp, ibi[ibi.size(0) - 1],
                           output->wav.hrv.aVLF);
  }
  // time values for interp.
  loop_ub = ibi.size(0);
  b_ibi.set_size(ibi.size(0));
  for (i = 0; i < loop_ub; i++) {
    b_ibi[i] = ibi[i];
  }
  loop_ub = ibi.size(0);
  c_ibi.set_size(ibi.size(0));
  for (i = 0; i < loop_ub; i++) {
    c_ibi[i] = ibi[i + ibi.size(0)];
  }
  coder::interp1(b_ibi, c_ibi, output->wav.hrv.aVLF, output->wav.hrv.LFHF);
  c_ibi.set_size(output->wav.hrv.LFHF.size(1));
  loop_ub = output->wav.hrv.LFHF.size(1);
  for (i = 0; i < loop_ub; i++) {
    c_ibi[i] = output->wav.hrv.LFHF[i];
  }
  // cubic spline interpolation
  // get limits of windows
  if (output->wav.hrv.aVLF[output->wav.hrv.aVLF.size(1) - 1] >= xtmp) {
    // windowArray - Seperates the input array s into smaller segments/windows
    // with lengths of winSize. winSize can be either in units of time or
    // samples. If winSize is given in units of time the first column of the
    // input array must containg time and method should be set to a value of 1.
    //
    // INPUTS:
    //    s: input array OR length of array. If outputing segs or using the time
    //        methods, the 1st dim must contain time values.
    //    winSize: size of segments to return
    //    method: method of segmenting array (0=by samples, 1=by time)
    //    overlap: amount that segments/windows overlap.
    //    flagCenters: decides whether to output center of windows instead of
    //    beginning
    //
    // OUTPUTS:
    //    idx: 2dim array that contains the start and end of each segment/window
    //    centers: array containg locations of center of windows (useful when
    //    using the time method)
    //
    // EXAMPLES:
    //    1. Seperate input using window size of 100 seconds and overlap
    //    of 50 seconds.
    //    idx=slidingWindow(input,100,50,1);
    //
    //    2. Seperate input using window size of 128 samples and overlap
    //    of 64 samples.
    //    idx=slidingWindow(input,128,64,0);
    //
    //    3. Seperate input using window size of 128 samples and overlap
    //    of 64 samples. Also outputs a cell array containing the actual
    //    segments [idx, segs]=slidingWindow(input,128,64,0,1);
    // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    //  Copyright (C) 2010, John T. Ramshur, jramshur@gmail.com
    //
    //  This file is part of HRVAS
    //
    //  HRVAS is free software: you can redistribute it and/or modify
    //  it under the terms of the GNU General Public License as published by
    //  the Free Software Foundation, either version 3 of the License, or
    //  (at your option) any later version.
    //
    //  HRVAS is distributed in the hope that it will be useful,
    //  but WITHOUT ANY WARRANTY; without even the implied warranty of
    //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    //  GNU General Public License for more details.
    //
    //  You should have received a copy of the GNU General Public License
    //  along with Foobar.  If not, see <http://www.gnu.org/licenses/>.
    // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    // Check inputs
    //  if nargin<5
    //      flagCenters=false;
    //      c=[];
    //      %includeSeg=0;
    //  elseif nargin<4
    //      overlap=0;
    //      includeSeg=0;
    //  elseif nargin<3
    //      overlap=0;
    //      method=0;
    //      includeSeg=0;
    //  elseif nargin<3
    //      error('slidingWindow: Too few input arguments.')
    //      return;
    //  end
    // check inputs
    //  if overlap>=winSize
    //      error('slidingWindow: overlap cannot be >= winSize.')
    //  end
    //  [m,n]=size(s);
    //  if (m==1) && (n==1) && (method==1)
    //      error('slidingWindow: Invalid input array for this windowing
    //      method.') return;
    //  end
    b_d = xtmp - overlap * fs;
    variance = (static_cast<double>(output->wav.hrv.aVLF.size(1)) - xtmp) + 1.0;
    if (std::isnan(b_d)) {
      iStart.set_size(1, 1);
      iStart[0] = rtNaN;
    } else if ((b_d == 0.0) || ((variance > 1.0) && (b_d < 0.0)) ||
               ((variance < 1.0) && (b_d > 0.0))) {
      iStart.set_size(1, 0);
    } else if (std::isinf(variance) && std::isinf(b_d)) {
      iStart.set_size(1, 1);
      iStart[0] = rtNaN;
    } else if (std::isinf(b_d)) {
      iStart.set_size(1, 1);
      iStart[0] = 1.0;
    } else if (std::floor(b_d) == b_d) {
      loop_ub = static_cast<int>((variance - 1.0) / b_d);
      iStart.set_size(1, loop_ub + 1);
      for (i = 0; i <= loop_ub; i++) {
        iStart[i] = b_d * static_cast<double>(i) + 1.0;
      }
    } else {
      coder::eml_float_colon(b_d, variance, iStart);
    }
    //      nseg = length(iStart);
    idx.set_size(iStart.size(1), 2);
    loop_ub = iStart.size(1);
    for (i = 0; i < loop_ub; i++) {
      idx[i] = iStart[i];
    }
    loop_ub = iStart.size(1);
    for (i = 0; i < loop_ub; i++) {
      idx[i + idx.size(0)] = iStart[i] + (xtmp - 1.0);
    }
    //  if includeSeg
    //      segs=cell(size(idx,1));
    //      for i=1:size(idx,1)
    //          segs{i}=s(idx(i,1):idx(i,2));%code here for outputing
    //      end
    //  end
    // (sample #)
  } else {
    idx.set_size(1, 2);
    idx[0] = 1.0;
    idx[1] = output->wav.hrv.aVLF.size(1);
  }
  xtmp = std::round(xtmp / 2.0);
  loop_ub = idx.size(0);
  output->ar.t.set_size(1, idx.size(0));
  for (i = 0; i < loop_ub; i++) {
    output->ar.t[i] = output->wav.hrv.aVLF[static_cast<int>(idx[i] + xtmp) - 1];
  }
  // calculate center time of window (s)
  // used for plotting
  // preallocate memory
  // number of PSD/windows
  nx = static_cast<int>(nfft);
  output->lomb.psd.set_size(nx, idx.size(0));
  output->ar.f.set_size(nx);
  for (i = 0; i < nx; i++) {
    output->ar.f[i] = 0.0;
  }
  // Calculate PSD
  i = idx.size(0);
  for (int b_i{0}; b_i < i; b_i++) {
    // Prepare y2 and t2
    d = idx[b_i];
    b_d = idx[b_i + idx.size(0)];
    if (d > b_d) {
      a1 = 0;
      y = -1;
    } else {
      a1 = static_cast<int>(d) - 1;
      y = static_cast<int>(b_d) - 1;
    }
    // remove linear trend
    //         y=detrend(y,'linear');
    loop_ub = y - a1;
    b_ibi.set_size(loop_ub + 1);
    for (y = 0; y <= loop_ub; y++) {
      b_ibi[y] = c_ibi[a1 + y];
    }
    b_d = coder::blockedSummation(b_ibi, loop_ub + 1) /
          (static_cast<double>(loop_ub) + 1.0);
    output->lomb.global.psd.set_size(loop_ub + 1);
    for (y = 0; y <= loop_ub; y++) {
      output->lomb.global.psd[y] = c_ibi[a1 + y] - b_d;
    }
    // remove mean
    loop_ub = output->lomb.global.psd.size(0);
    coder::hamming(static_cast<double>(output->lomb.global.psd.size(0)), psd);
    if (output->lomb.global.psd.size(0) == psd.size(0)) {
      for (a1 = 0; a1 < loop_ub; a1++) {
        output->lomb.global.psd[a1] = output->lomb.global.psd[a1] * psd[a1];
      }
    } else {
      binary_expand_op(output, psd);
    }
    // hamming window
    // Calculate PSD
    coder::pburg(output->lomb.global.psd, AR_order, nfft * 2.0 - 1.0, fs, psd,
                 output->ar.f);
    loop_ub = psd.size(0);
    for (a1 = 0; a1 < loop_ub; a1++) {
      output->lomb.psd[a1 + output->lomb.psd.size(0) * b_i] = psd[a1];
    }
  }
  output->ar.psd.set_size(output->lomb.psd.size(0), output->lomb.psd.size(1));
  loop_ub = output->lomb.psd.size(0) * output->lomb.psd.size(1);
  for (i = 0; i < loop_ub; i++) {
    output->ar.psd[i] = output->lomb.psd[i];
  }
  calcHRV(output->ar.f, output->lomb.psd, VLF, LF, HF, &output->ar.hrv);
  // global psd
  output->ar.global.f.set_size(output->ar.f.size(0));
  loop_ub = output->ar.f.size(0);
  for (i = 0; i < loop_ub; i++) {
    output->ar.global.f[i] = output->ar.f[i];
  }
  coder::mean(output->lomb.psd, output->lomb.global.psd);
  output->ar.global.psd.set_size(output->lomb.global.psd.size(0));
  loop_ub = output->lomb.global.psd.size(0);
  for (i = 0; i < loop_ub; i++) {
    output->ar.global.psd[i] = output->lomb.global.psd[i];
  }
  calcAreas(output->ar.f, output->lomb.global.psd, VLF, LF, HF,
            &output->ar.global.hrv);
  // Lomb
  //      output.lomb.t=t;
  loop_ub = ibi.size(0);
  b_ibi.set_size(ibi.size(0));
  for (i = 0; i < loop_ub; i++) {
    b_ibi[i] = ibi[i];
  }
  loop_ub = ibi.size(0);
  c_ibi.set_size(ibi.size(0));
  for (i = 0; i < loop_ub; i++) {
    c_ibi[i] = ibi[i + ibi.size(0)];
  }
  calcLomb(b_ibi, c_ibi, nfft, fs / 2.0, winSize, overlap, output->lomb.psd,
           psd, output->lomb.t);
  output->lomb.f.set_size(psd.size(0));
  loop_ub = psd.size(0);
  for (i = 0; i < loop_ub; i++) {
    output->lomb.f[i] = psd[i];
  }
  calcHRV(psd, output->lomb.psd, VLF, LF, HF, &output->lomb.hrv);
  // global psd
  output->lomb.global.f.set_size(psd.size(0));
  loop_ub = psd.size(0);
  for (i = 0; i < loop_ub; i++) {
    output->lomb.global.f[i] = psd[i];
  }
  coder::mean(output->lomb.psd, output->lomb.global.psd);
  calcAreas(psd, output->lomb.global.psd, VLF, LF, HF,
            &output->lomb.global.hrv);
  // Wavelet
  // y=nibi(:,2);
  //          clear nibi; % don't need it anymore
  if (std::isnan(ibi[0]) || std::isnan(d_tmp) ||
      std::isnan(ibi[ibi.size(0) - 1])) {
    iStart.set_size(1, 1);
    iStart[0] = rtNaN;
  } else if ((d_tmp == 0.0) ||
             ((ibi[0] < ibi[ibi.size(0) - 1]) && (d_tmp < 0.0)) ||
             ((ibi[ibi.size(0) - 1] < ibi[0]) && (d_tmp > 0.0))) {
    iStart.set_size(1, 0);
  } else if ((std::isinf(ibi[0]) || std::isinf(ibi[ibi.size(0) - 1])) &&
             (std::isinf(d_tmp) || (ibi[0] == ibi[ibi.size(0) - 1]))) {
    iStart.set_size(1, 1);
    iStart[0] = rtNaN;
  } else if (std::isinf(d_tmp)) {
    iStart.set_size(1, 1);
    iStart[0] = ibi[0];
  } else if ((std::floor(ibi[0]) == ibi[0]) && (std::floor(d_tmp) == d_tmp)) {
    iStart.set_size(
        1, static_cast<int>((ibi[ibi.size(0) - 1] - ibi[0]) / d_tmp) + 1);
    loop_ub = static_cast<int>((ibi[ibi.size(0) - 1] - ibi[0]) / d_tmp);
    for (i = 0; i <= loop_ub; i++) {
      iStart[i] = ibi[0] + d_tmp * static_cast<double>(i);
    }
  } else {
    coder::eml_float_colon(ibi[0], d_tmp, ibi[ibi.size(0) - 1], iStart);
  }
  // time values for interp.
  loop_ub = ibi.size(0);
  b_ibi.set_size(ibi.size(0));
  for (i = 0; i < loop_ub; i++) {
    b_ibi[i] = ibi[i];
  }
  loop_ub = ibi.size(0);
  c_ibi.set_size(ibi.size(0));
  for (i = 0; i < loop_ub; i++) {
    c_ibi[i] = ibi[i + ibi.size(0)];
  }
  coder::interp1(b_ibi, c_ibi, iStart, output->wav.hrv.LFHF);
  c_ibi.set_size(output->wav.hrv.LFHF.size(1));
  loop_ub = output->wav.hrv.LFHF.size(1);
  for (i = 0; i < loop_ub; i++) {
    c_ibi[i] = output->wav.hrv.LFHF[i];
  }
  // cubic spline interpolation
  output->wav.t.set_size(1, iStart.size(1));
  loop_ub = iStart.size(1);
  for (i = 0; i < loop_ub; i++) {
    output->wav.t[i] = iStart[i];
  }
  xtmp = coder::b_std(c_ibi);
  xtmp = std::sqrt(xtmp * xtmp);
  b_d = coder::blockedSummation(c_ibi, c_ibi.size(0)) /
        static_cast<double>(c_ibi.size(0));
  psd.set_size(c_ibi.size(0));
  loop_ub = c_ibi.size(0);
  for (i = 0; i < loop_ub; i++) {
    psd[i] = (c_ibi[i] - b_d) / xtmp;
  }
  // xlim = [0,t2(end)];  % plotting range
  //  pad the time series with zeroes (recommended)
  //  this will do 4 sub-octaves per octave
  xtmp = 2.0 * d_tmp;
  //  this says start at a scale of 0.5 s
  //  this says do 7 powers-of-two with dj sub-octaves each
  //  lag-1 autocorrelation for red noise background
  // mother = 'DOG';
  // mother = 'Paul';
  //  this is for the MORLET wavelet
  //  Wavelet transform
  // WAVELET  1D Wavelet transform with optional singificance testing
  //
  //    [WAVE,PERIOD,SCALE,COI] = wavelet(Y,DT,PAD,DJ,S0,J1,MOTHER,PARAM)
  //
  //    Computes the wavelet transform of the vector Y (length N),
  //    with sampling rate DT.
  //
  //    By default, the Morlet wavelet (k0=6) is used.
  //    The wavelet basis is normalized to have total energy=1 at all scales.
  //
  //  INPUTS:
  //
  //     Y = the time series of length N.
  //     DT = amount of time between each Y value, i.e. the sampling time.
  //
  //  OUTPUTS:
  //
  //     WAVE is the WAVELET transform of Y. This is a complex array
  //     of dimensions (N,J1+1). FLOAT(WAVE) gives the WAVELET amplitude,
  //     ATAN(IMAGINARY(WAVE),FLOAT(WAVE) gives the WAVELET phase.
  //     The WAVELET power spectrum is ABS(WAVE)^2.
  //     Its units are sigma^2 (the time series variance).
  //
  //  OPTIONAL INPUTS:
  //
  //  *** Note *** setting any of the following to -1 will cause the default
  //                value to be used.
  //
  //     PAD = if set to 1 (default is 0), pad time series with enough zeroes to
  //     get
  //          N up to the next higher power of 2. This prevents wraparound
  //          from the end of the time series to the beginning, and also
  //          speeds up the FFT's used to do the wavelet transform.
  //          This will not eliminate all edge effects (see COI below).
  //
  //     DJ = the spacing between discrete scales. Default is 0.25.
  //          A smaller # will give better scale resolution, but be slower to
  //          plot.
  //
  //     S0 = the smallest scale of the wavelet.  Default is 2*DT.
  //
  //     J1 = the # of scales minus one. Scales range from S0 up to
  //     S0*2^(J1*DJ),
  //         to give a total of (J1+1) scales. Default is J1 = (LOG2(N
  //         DT/S0))/DJ.
  //
  //     MOTHER = the mother wavelet function.
  //              The choices are 'MORLET', 'PAUL', or 'DOG'
  //
  //     PARAM = the mother wavelet parameter.
  //             For 'MORLET' this is k0 (wavenumber), default is 6.
  //             For 'PAUL' this is m (order), default is 4.
  //             For 'DOG' this is m (m-th derivative), default is 2.
  //
  //  OPTIONAL OUTPUTS:
  //
  //     PERIOD = the vector of "Fourier" periods (in time units) that
  //     corresponds
  //            to the SCALEs.
  //
  //     SCALE = the vector of scale indices, given by S0*2^(j*DJ), j=0...J1
  //             where J1+1 is the total # of scales.
  //
  //     COI = if specified, then return the Cone-of-Influence, which is a
  //     vector
  //         of N points that contains the maximum period of useful information
  //         at that particular time.
  //         Periods greater than this are subject to edge effects.
  //         This can be used to plot COI lines on a contour plot by doing:
  //
  //               contour(time,log(period),log(power))
  //               plot(time,log(coi),'k')
  //
  // ----------------------------------------------------------------------------
  //    Copyright (C) 1995-2004, Christopher Torrence and Gilbert P. Compo
  //
  //    This software may be used, copied, or redistributed as long as it is not
  //    sold and this copyright notice is reproduced on each copy made. This
  //    routine is provided as is without any express or implied warranties
  //    whatsoever.
  //
  //  Notice: Please acknowledge the use of the above software in any
  //  publications:
  //     ``Wavelet software was provided by C. Torrence and G. Compo,
  //       and is available at URL:
  //       http://paos.colorado.edu/research/wavelets/''.
  //
  //  Reference: Torrence, C. and G. P. Compo, 1998: A Practical Guide to
  //             Wavelet Analysis. <I>Bull. Amer. Meteor. Soc.</I>, 79, 61-78.
  //
  //  Please send a copy of such publications to either C. Torrence or G. Compo:
  //   Dr. Christopher Torrence               Dr. Gilbert P. Compo
  //   Research Systems, Inc.                 Climate Diagnostics Center
  //   4990 Pearl East Circle                 325 Broadway R/CDC1
  //   Boulder, CO 80301, USA                 Boulder, CO 80305-3328, USA
  //   E-mail: chris[AT]rsinc[DOT]com         E-mail: compo[AT]colorado[DOT]edu
  // ----------------------------------------------------------------------------
  // ....construct time series to analyze, pad if necessary
  iStart.set_size(1, psd.size(0));
  loop_ub = psd.size(0);
  for (i = 0; i < loop_ub; i++) {
    iStart[i] = 0.0;
  }
  b_d = coder::blockedSummation(psd, psd.size(0)) /
        static_cast<double>(psd.size(0));
  if (psd.size(0) < 1) {
    loop_ub = 0;
  } else {
    loop_ub = psd.size(0);
  }
  for (i = 0; i < loop_ub; i++) {
    iStart[i] = psd[i] - b_d;
  }
  //  power of 2 nearest to N
  i = iStart.size(1);
  loop_ub = static_cast<int>(
      rt_powd_snf(2.0, std::trunc(std::log(static_cast<double>(psd.size(0))) /
                                      0.69314718055994529 +
                                  0.4999) +
                           1.0) -
      static_cast<double>(psd.size(0)));
  iStart.set_size(
      iStart.size(0),
      iStart.size(1) +
          static_cast<int>(
              rt_powd_snf(
                  2.0, std::trunc(std::log(static_cast<double>(psd.size(0))) /
                                      0.69314718055994529 +
                                  0.4999) +
                           1.0) -
              static_cast<double>(psd.size(0))));
  for (a1 = 0; a1 < loop_ub; a1++) {
    iStart[i + a1] = 0.0;
  }
  // ....construct wavenumber array used in transform [Eqn(5)]
  nx = static_cast<int>(std::trunc(static_cast<double>(iStart.size(1)) / 2.0));
  if (nx < 1) {
    output->wav.hrv.aVLF.set_size(output->wav.hrv.aVLF.size(0), 0);
  } else {
    output->wav.hrv.aVLF.set_size(1, nx);
    loop_ub = nx - 1;
    for (i = 0; i <= loop_ub; i++) {
      output->wav.hrv.aVLF[i] = static_cast<double>(i) + 1.0;
    }
  }
  output->wav.hrv.aVLF.set_size(1, output->wav.hrv.aVLF.size(1));
  d = 6.2831853071795862 / (static_cast<double>(iStart.size(1)) * d_tmp);
  loop_ub = output->wav.hrv.aVLF.size(1) - 1;
  for (i = 0; i <= loop_ub; i++) {
    output->wav.hrv.aVLF[i] = output->wav.hrv.aVLF[i] * d;
  }
  i = static_cast<int>(
      std::trunc((static_cast<double>(iStart.size(1)) - 1.0) / 2.0));
  if (i < 1) {
    i = 0;
    a1 = 1;
    y = -1;
  } else {
    i--;
    a1 = -1;
    y = 0;
  }
  loop_ub = div_s32(y - i, a1);
  r.set_size(1, (output->wav.hrv.aVLF.size(1) + loop_ub) + 2);
  r[0] = 0.0;
  nx = output->wav.hrv.aVLF.size(1);
  for (y = 0; y < nx; y++) {
    r[y + 1] = output->wav.hrv.aVLF[y];
  }
  for (y = 0; y <= loop_ub; y++) {
    r[(y + output->wav.hrv.aVLF.size(1)) + 1] =
        -output->wav.hrv.aVLF[i + a1 * y];
  }
  output->wav.hrv.aVLF.set_size(1, r.size(1));
  loop_ub = r.size(1);
  for (i = 0; i < loop_ub; i++) {
    output->wav.hrv.aVLF[i] = r[i];
  }
  // ....compute FFT of the (padded) time series
  coder::fft(iStart, f);
  //  [Eqn(3)]
  // ....construct SCALE array & empty PERIOD & WAVE arrays
  for (i = 0; i < 449; i++) {
    output->wav.global.psd[i] = xtmp * b[i];
  }
  //  define the wavelet array
  wave.set_size(449, iStart.size(1));
  loop_ub = 449 * iStart.size(1);
  for (i = 0; i < loop_ub; i++) {
    wave[i].re = 0.0;
    wave[i].im = 0.0;
  }
  //  make it complex
  //  loop through all scales and compute transform
  loop_ub = output->wav.hrv.aVLF.size(1);
  for (a1 = 0; a1 < 449; a1++) {
    // WAVE_BASES  1D Wavelet functions Morlet, Paul, or DOG
    //
    //   [DAUGHTER,FOURIER_FACTOR,COI,DOFMIN] = ...
    //       wave_bases(MOTHER,K,SCALE,PARAM);
    //
    //    Computes the wavelet function as a function of Fourier frequency,
    //    used for the wavelet transform in Fourier space.
    //    (This program is called automatically by WAVELET)
    //
    //  INPUTS:
    //
    //     MOTHER = a string, equal to 'MORLET' or 'PAUL' or 'DOG'
    //     K = a vector, the Fourier frequencies at which to calculate the
    //     wavelet SCALE = a number, the wavelet scale PARAM = the
    //     nondimensional parameter for the wavelet function
    //
    //  OUTPUTS:
    //
    //     DAUGHTER = a vector, the wavelet function
    //     FOURIER_FACTOR = the ratio of Fourier period to scale
    //     COI = a number, the cone-of-influence size at the scale
    //     DOFMIN = a number, degrees of freedom for each point in the wavelet
    //     power
    //              (either 2 for Morlet and Paul, or 1 for the DOG)
    //
    // ----------------------------------------------------------------------------
    //    Copyright (C) 1995-1998, Christopher Torrence and Gilbert P. Compo
    //    University of Colorado, Program in Atmospheric and Oceanic Sciences.
    //    This software may be used, copied, or redistributed as long as it is
    //    not sold and this copyright notice is reproduced on each copy made.
    //    This routine is provided as is without any express or implied
    //    warranties whatsoever.
    // ----------------------------------------------------------------------------
    // -----------------------------------  Morlet
    iStart.set_size(1, output->wav.hrv.aVLF.size(1));
    for (i = 0; i < loop_ub; i++) {
      xtmp = output->wav.global.psd[a1] * output->wav.hrv.aVLF[i] - 6.0;
      iStart[i] = xtmp * xtmp;
    }
    if (iStart.size(1) == output->wav.hrv.aVLF.size(1)) {
      nx = iStart.size(1) - 1;
      iStart.set_size(1, iStart.size(1));
      for (i = 0; i <= nx; i++) {
        iStart[i] = -iStart[i] / 2.0 *
                    static_cast<double>(output->wav.hrv.aVLF[i] > 0.0);
      }
    } else {
      binary_expand_op(iStart, output);
    }
    xtmp = std::sqrt(output->wav.global.psd[a1] * output->wav.hrv.aVLF[1]) *
           0.75112554446494251 *
           std::sqrt(static_cast<double>(output->wav.hrv.aVLF.size(1)));
    //  total energy=N   [Eqn(7)]
    nx = iStart.size(1);
    for (int b_i{0}; b_i < nx; b_i++) {
      iStart[b_i] = std::exp(iStart[b_i]);
    }
    //  Heaviside step function
    //  Scale-->Fourier [Sec.3h]
    //  Cone-of-influence [Sec.3g]
    //  Degrees of freedom
    if (iStart.size(1) == 1) {
      i = output->wav.hrv.aVLF.size(1);
    } else {
      i = iStart.size(1);
    }
    if ((iStart.size(1) == output->wav.hrv.aVLF.size(1)) && (f.size(1) == i)) {
      b_f.set_size(1, f.size(1));
      nx = f.size(1);
      for (i = 0; i < nx; i++) {
        d = xtmp * iStart[i];
        b_f[i].re =
            d * static_cast<double>(output->wav.hrv.aVLF[i] > 0.0) * f[i].re;
        b_f[i].im =
            d * static_cast<double>(output->wav.hrv.aVLF[i] > 0.0) * f[i].im;
      }
      coder::ifft(b_f, r2);
      nx = r2.size(1);
      for (i = 0; i < nx; i++) {
        wave[a1 + 449 * i] = r2[i];
      }
    } else {
      binary_expand_op(wave, a1, f, xtmp, iStart, output);
    }
    //  wavelet transform[Eqn(4)]
  }
  //  COI [Sec.3g]
  if (psd.size(0) < 1) {
    i = 0;
  } else {
    i = psd.size(0);
  }
  if (psd.size(0) < 1) {
    loop_ub = 0;
  } else {
    loop_ub = psd.size(0);
  }
  wave.set_size(449, loop_ub);
  //  get rid of padding before returning
  //  Reference: Torrence, C. and G. P. Compo, 1998: A Practical Guide to
  //  Wavelet Analysis. <I>Bull. Amer. Meteor. Soc.</I>, 79, 61-78.
  nx = 449 * i;
  output->wav.psd.set_size(449, i);
  for (int b_i{0}; b_i < nx; b_i++) {
    output->wav.psd[b_i] = rt_hypotd_snf(wave[b_i].re, wave[b_i].im);
  }
  loop_ub = 449 * output->wav.psd.size(1);
  output->wav.psd.set_size(449, output->wav.psd.size(1));
  for (i = 0; i < loop_ub; i++) {
    xtmp = output->wav.psd[i];
    output->wav.psd[i] = rt_powd_snf(xtmp, 2.0);
  }
  //  compute wavelet power spectrum
  for (i = 0; i < 449; i++) {
    output->wav.f[i] = 1.0 / (1.0330436477492537 * output->wav.global.psd[i]);
  }
  for (nx = 0; nx < 224; nx++) {
    xtmp = output->wav.f[nx];
    output->wav.f[nx] = output->wav.f[448 - nx];
    output->wav.f[448 - nx] = xtmp;
  }
  // frequency in ascending order
  nx = output->wav.psd.size(1);
  for (a1 = 0; a1 < nx; a1++) {
    for (int b_i{0}; b_i < 224; b_i++) {
      xtmp = output->wav.psd[b_i + 449 * a1];
      output->wav.psd[b_i + 449 * a1] = output->wav.psd[(449 * a1 - b_i) + 448];
      output->wav.psd[(449 * a1 - b_i) + 448] = xtmp;
    }
  }
  // flip to match freq. order
  variance = coder::var(c_ibi);
  // calcAreas - Calulates areas/energy under the PSD curve within the freq
  // bands defined by VLF, LF, and HF. Returns areas/energies as ms^2,
  // percentage, and normalized units. Also returns LF/HF ratio.
  //
  // Inputs:
  //    PSD: PSD vector
  //    F: Freq vector
  //    VLF, LF, HF: array containing VLF, LF, and HF freq limits
  //    flagNormalize: option to normalize PSD to max(PSD)
  // Output:
  //
  // Usage:
  //
  //  Scale-average between VLF, LF, and HF bands/scales
  // f=fliplr(f); %put f in it's original order
  //  expand scale --> (J+1)x(N) array
  r1.set_size(449, c_ibi.size(0));
  loop_ub = c_ibi.size(0);
  for (i = 0; i < loop_ub; i++) {
    for (a1 = 0; a1 < 449; a1++) {
      r1[a1 + 449 * i] = output->wav.global.psd[a1];
    }
  }
  if (output->wav.psd.size(1) == r1.size(1)) {
    scale_avg.set_size(449, c_ibi.size(0));
    loop_ub = c_ibi.size(0);
    for (i = 0; i < loop_ub; i++) {
      for (a1 = 0; a1 < 449; a1++) {
        scale_avg[a1 + 449 * i] =
            output->wav.psd[a1 + 449 * i] / output->wav.global.psd[a1];
      }
    }
  } else {
    binary_expand_op(scale_avg, output, c_ibi);
  }
  //  [Eqn(24)]
  xtmp = variance * 0.015625 * d_tmp / 0.776;
  d = VLF[0];
  b_d = VLF[1];
  for (i = 0; i < 449; i++) {
    b_output[i] = ((output->wav.f[i] >= d) && (output->wav.f[i] < b_d));
  }
  coder::eml_find(b_output, tmp_data, tmp_size);
  loop_ub = scale_avg.size(1);
  b_scale_avg.set_size(tmp_size[1], scale_avg.size(1));
  for (i = 0; i < loop_ub; i++) {
    nx = tmp_size[1];
    for (a1 = 0; a1 < nx; a1++) {
      b_scale_avg[a1 + b_scale_avg.size(0) * i] =
          scale_avg[(tmp_data[a1] + 449 * i) - 1];
    }
  }
  coder::sum(b_scale_avg, output->wav.hrv.aVLF);
  output->wav.hrv.aVLF.set_size(1, output->wav.hrv.aVLF.size(1));
  loop_ub = output->wav.hrv.aVLF.size(1) - 1;
  for (i = 0; i <= loop_ub; i++) {
    output->wav.hrv.aVLF[i] = xtmp * output->wav.hrv.aVLF[i];
  }
  //  [Eqn(24)]
  d = LF[0];
  b_d = LF[1];
  for (i = 0; i < 449; i++) {
    b_output[i] = ((output->wav.f[i] >= d) && (output->wav.f[i] < b_d));
  }
  coder::eml_find(b_output, tmp_data, tmp_size);
  loop_ub = scale_avg.size(1);
  b_scale_avg.set_size(tmp_size[1], scale_avg.size(1));
  for (i = 0; i < loop_ub; i++) {
    nx = tmp_size[1];
    for (a1 = 0; a1 < nx; a1++) {
      b_scale_avg[a1 + b_scale_avg.size(0) * i] =
          scale_avg[(tmp_data[a1] + 449 * i) - 1];
    }
  }
  coder::sum(b_scale_avg, output->wav.hrv.LFHF);
  output->wav.hrv.LFHF.set_size(1, output->wav.hrv.LFHF.size(1));
  loop_ub = output->wav.hrv.LFHF.size(1) - 1;
  for (i = 0; i <= loop_ub; i++) {
    output->wav.hrv.LFHF[i] = xtmp * output->wav.hrv.LFHF[i];
  }
  d = HF[0];
  b_d = HF[1];
  for (i = 0; i < 449; i++) {
    b_output[i] = ((output->wav.f[i] >= d) && (output->wav.f[i] < b_d));
  }
  coder::eml_find(b_output, tmp_data, tmp_size);
  loop_ub = scale_avg.size(1);
  b_scale_avg.set_size(tmp_size[1], scale_avg.size(1));
  for (i = 0; i < loop_ub; i++) {
    nx = tmp_size[1];
    for (a1 = 0; a1 < nx; a1++) {
      b_scale_avg[a1 + b_scale_avg.size(0) * i] =
          scale_avg[(tmp_data[a1] + 449 * i) - 1];
    }
  }
  coder::sum(b_scale_avg, output->wav.hrv.aHF);
  output->wav.hrv.aHF.set_size(1, output->wav.hrv.aHF.size(1));
  loop_ub = output->wav.hrv.aHF.size(1) - 1;
  for (i = 0; i <= loop_ub; i++) {
    output->wav.hrv.aHF[i] = xtmp * output->wav.hrv.aHF[i];
  }
  //  calculate raw areas (power under curve), within the freq bands (ms^2)
  output->wav.hrv.aLF.set_size(1, output->wav.hrv.LFHF.size(1));
  loop_ub = output->wav.hrv.LFHF.size(1) - 1;
  for (i = 0; i <= loop_ub; i++) {
    output->wav.hrv.aLF[i] = output->wav.hrv.LFHF[i];
  }
  if (output->wav.hrv.aVLF.size(1) == 1) {
    i = output->wav.hrv.LFHF.size(1);
  } else {
    i = output->wav.hrv.aVLF.size(1);
  }
  if ((output->wav.hrv.aVLF.size(1) == output->wav.hrv.LFHF.size(1)) &&
      (i == output->wav.hrv.aHF.size(1))) {
    iStart.set_size(1, output->wav.hrv.aVLF.size(1));
    loop_ub = output->wav.hrv.aVLF.size(1);
    for (i = 0; i < loop_ub; i++) {
      iStart[i] = (output->wav.hrv.aVLF[i] + output->wav.hrv.LFHF[i]) +
                  output->wav.hrv.aHF[i];
    }
  } else {
    c_binary_expand_op(iStart, output);
  }
  output->wav.hrv.aTotal.set_size(1, iStart.size(1));
  loop_ub = iStart.size(1);
  for (i = 0; i < loop_ub; i++) {
    output->wav.hrv.aTotal[i] = iStart[i];
  }
  // calculate areas relative to the total area (%)
  if (output->wav.hrv.aVLF.size(1) == iStart.size(1)) {
    output->wav.hrv.pVLF.set_size(1, output->wav.hrv.aVLF.size(1));
    loop_ub = output->wav.hrv.aVLF.size(1);
    for (i = 0; i < loop_ub; i++) {
      output->wav.hrv.pVLF[i] = output->wav.hrv.aVLF[i] / iStart[i] * 100.0;
    }
  } else {
    e_binary_expand_op(output, iStart);
  }
  if (output->wav.hrv.LFHF.size(1) == iStart.size(1)) {
    output->wav.hrv.pLF.set_size(1, output->wav.hrv.LFHF.size(1));
    loop_ub = output->wav.hrv.LFHF.size(1);
    for (i = 0; i < loop_ub; i++) {
      output->wav.hrv.pLF[i] = output->wav.hrv.LFHF[i] / iStart[i] * 100.0;
    }
  } else {
    d_binary_expand_op(output, iStart);
  }
  if (output->wav.hrv.aHF.size(1) == iStart.size(1)) {
    output->wav.hrv.pHF.set_size(1, output->wav.hrv.aHF.size(1));
    loop_ub = output->wav.hrv.aHF.size(1);
    for (i = 0; i < loop_ub; i++) {
      output->wav.hrv.pHF[i] = output->wav.hrv.aHF[i] / iStart[i] * 100.0;
    }
  } else {
    c_binary_expand_op(output, iStart);
  }
  // calculate normalized areas (relative to HF+LF, n.u.)
  if (output->wav.hrv.LFHF.size(1) == output->wav.hrv.aHF.size(1)) {
    iStart.set_size(1, output->wav.hrv.LFHF.size(1));
    loop_ub = output->wav.hrv.LFHF.size(1);
    for (i = 0; i < loop_ub; i++) {
      iStart[i] = output->wav.hrv.LFHF[i] + output->wav.hrv.aHF[i];
    }
  } else {
    b_binary_expand_op(iStart, output);
  }
  if (output->wav.hrv.LFHF.size(1) == iStart.size(1)) {
    output->wav.hrv.nLF.set_size(1, output->wav.hrv.LFHF.size(1));
    loop_ub = output->wav.hrv.LFHF.size(1);
    for (i = 0; i < loop_ub; i++) {
      output->wav.hrv.nLF[i] = output->wav.hrv.LFHF[i] / iStart[i];
    }
  } else {
    b_binary_expand_op(output, iStart);
  }
  if (output->wav.hrv.aHF.size(1) == iStart.size(1)) {
    output->wav.hrv.nHF.set_size(1, output->wav.hrv.aHF.size(1));
    loop_ub = output->wav.hrv.aHF.size(1);
    for (i = 0; i < loop_ub; i++) {
      output->wav.hrv.nHF[i] = output->wav.hrv.aHF[i] / iStart[i];
    }
  } else {
    binary_expand_op(output, iStart);
  }
  // calculate LF/HF ratio
  if (output->wav.hrv.LFHF.size(1) == output->wav.hrv.aHF.size(1)) {
    loop_ub = output->wav.hrv.LFHF.size(1) - 1;
    output->wav.hrv.LFHF.set_size(1, output->wav.hrv.LFHF.size(1));
    for (i = 0; i <= loop_ub; i++) {
      output->wav.hrv.LFHF[i] =
          output->wav.hrv.LFHF[i] / output->wav.hrv.aHF[i];
    }
  } else {
    binary_expand_op(output);
  }
  x.set_size(1, output->wav.hrv.LFHF.size(1));
  loop_ub = output->wav.hrv.LFHF.size(1);
  for (i = 0; i < loop_ub; i++) {
    x[i] = (output->wav.hrv.LFHF[i] > 1.0);
  }
  nx = x.size(1);
  if (x.size(1) == 0) {
    y = 0;
  } else {
    y = x[0];
    for (int b_i{2}; b_i <= nx; b_i++) {
      y += x[b_i - 1];
    }
  }
  x.set_size(1, output->wav.hrv.LFHF.size(1));
  loop_ub = output->wav.hrv.LFHF.size(1);
  for (i = 0; i < loop_ub; i++) {
    x[i] = (output->wav.hrv.LFHF[i] <= 1.0);
  }
  nx = x.size(1);
  if (x.size(1) == 0) {
    a1 = 0;
  } else {
    a1 = x[0];
    for (int b_i{2}; b_i <= nx; b_i++) {
      a1 += x[b_i - 1];
    }
  }
  output->wav.hrv.rLFHF = static_cast<double>(y) / static_cast<double>(a1);
  // calculate peaks
  output->wav.hrv.peakVLF.set_size(1, output->wav.hrv.aVLF.size(1));
  loop_ub = output->wav.hrv.aVLF.size(1);
  for (i = 0; i < loop_ub; i++) {
    output->wav.hrv.peakVLF[i] = 0.0;
  }
  output->wav.hrv.peakLF.set_size(1, output->wav.hrv.aVLF.size(1));
  loop_ub = output->wav.hrv.aVLF.size(1);
  for (i = 0; i < loop_ub; i++) {
    output->wav.hrv.peakLF[i] = 0.0;
  }
  output->wav.hrv.peakHF.set_size(1, output->wav.hrv.aVLF.size(1));
  loop_ub = output->wav.hrv.aVLF.size(1);
  for (i = 0; i < loop_ub; i++) {
    output->wav.hrv.peakHF[i] = 0.0;
  }
  //  Global wavelet power spectrum
  coder::sum(output->wav.psd, output->wav.global.psd);
  for (i = 0; i < 449; i++) {
    output->wav.global.psd[i] = variance * (output->wav.global.psd[i] /
                                            static_cast<double>(c_ibi.size(0)));
    output->wav.global.f[i] = output->wav.f[i];
  }
  calcAreas(output->wav.f, output->wav.global.psd, VLF, LF, HF,
            &output->wav.global.hrv);
}

//
// File trailer for timeFreqHRV.cpp
//
// [EOF]
//
